require(Rcpp)
output_version <- '2020-10-20'

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01') %>%
  st_set_crs(.,4326)
ntl.extent <- extent(as(ntl.zones,"Spatial"))

ntl.zones.prj <- ntl.zones %>%
  # 102022
  st_transform(crs="+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs") %>% 
  st_centroid() %>% 
  st_geometry()

M.ntl.zones <- ntl.zones.prj %>% 
  st_transform(4326) %>%
  st_coordinates()

# convert polygons to points 
ntl.zones.pts = st_as_sf(as.data.frame(M.ntl.zones), coords = c("X", "Y"), crs = 4326)
ntl.zones = st_drop_geometry(subset(ntl.zones,select=c(OBJECTID, GID_0)))
ntl.zones.pts = dplyr::bind_cols(ntl.zones.pts,ntl.zones)

ntl.zones <- read_sf(dsn=paste0(wkdir,'/local/gis_data/003_boundaries/fishnet'),layer='fishnet_lake_chad_ntl_adm0_x_d01')

afclima_fn <- file.path(afclima_dir,'ClimAf_1_1981_2050_Suitweight_HIS95_12crops.tif')
afclima <- raster(afclima_fn)

# manually set outside for 5x5
ntl.extent <- c(-10,35,-10,35)
afclima.crop <- crop(afclima,ntl.extent,snap='out')

## - create measure avg same cell, avg 5 cells from cell, avg 10 cells from cell
afclima.diam.crop <-  sqrt(raster::area(afclima.crop))
decdeg050km = 50 / cellStats(afclima.diam.crop,'mean')
lendd050km = res(afclima.crop)[1] * decdeg050km

decdeg100km = 100 / cellStats(afclima.diam.crop,'mean')
lendd100km = res(afclima.crop)[1] * decdeg100km

fw5x5 <- focalWeight(afclima.crop, lendd050km, type='circle')
fw10x10 <- focalWeight(afclima.crop, lendd100km, type='circle')

agsuitclrad050 <- focal(afclima.crop, w=fw5x5, fun=sum,na.rm=TRUE)
agsuitclrad100 <- focal(afclima.crop, w=fw10x10, fun=sum,na.rm=TRUE)

ntl.zones$agsuit <- exactextractr::exact_extract(afclima.crop,ntl.zones,fun='mean')
ntl.zones$agsuit050 <- exactextractr::exact_extract(agsuitclrad050,ntl.zones,fun='mean')
ntl.zones$agsuit100 <- exactextractr::exact_extract(agsuitclrad100,ntl.zones,fun='mean')

save.dta13(as.data.frame(st_drop_geometry(subset(ntl.zones,select=c(OBJECTID,GID_0,agsuit,agsuit050,agsuit100)))), 
           paste0(wkdir,"/proc_data/AGSUIT",output_version,".dta"))